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Abstract 

We consider suspensions of neutrally-buoyant finite-size rigid spherical particles 
in channel flo'w and investigate the relation between the particle dynamics and 
the mean bulk behavior of the mixture for Reynolds numbers 500 < Re < 5000 
and particle volume fraction 0 < <f> < 0.3, via fully resolved numerical sim¬ 
ulations. Analysis of the momentum balance reveals the existence of three 
different regimes: laminar, turbulent and inertial shear-thickening depending 
on which of the stress terms, viscous, Reynolds or particle stress, is the major 
responsible for the momentum transfer across the channel. We show that both 
Reynolds and particle stress dominated flows fall into the Bagnoldian inertial 
regime and that the Bagnold number can predict the bulk behavior although 
this is due to two distinct physical mechanisms. A turbulent flow is charac¬ 
terized by larger particle dispersion and a more uniform particle distribution, 
whereas the particulate-dominated flows is associated with a significant particle 
migration towards the channel center where the flow is smooth laminar-like and 
dispersion low. Interestingly, the collision kernel shows similar values in the 
different regimes, although the relative particle velocity and clustering clearly 
vary with inertia and particle concentration. 
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1. Introduction 


Particles suspended in a carrier fluid can be found in many biological, geo¬ 
physical and industrial flows. Some obvious examples are the blood flow in the 
human body, pyroclastic flows from volcanos, sedimentations in sea beds, flu¬ 
idized beds and slurry flows. Moreover, the knowledge of the particle dynamics 
is relevant, among others, in biomechanical applications for extracorporeal de¬ 
vices and formation of clots. Suspensions are typically employed to transport 
and mix particles by means of a carrier fluid [li. The overall effect of particles 
on the flow dynamics has therefore a significant impact on the energy consump¬ 
tion of biological and industrial processes. Despite the numerous applications, 
however, it is still difficult to estimate the force needed to drive suspensions 
and the internal dissipation mechanisms are not fully understood, especially 
in a turbulent flow. Unlike single phase flows where the pressure drop can be 
accurately predicted as a function of the Reynolds number and the properties 
of the wall surface (roughness effects), additional parameters become relevant 
in the presence of a suspended phase when the properties of the particles (size, 
shape, density, stiffness, volume fraction, mass fraction) affect the overall dy¬ 
namics of the suspension. The behavior of these multiphase flows becomes even 
more complicated when the particle volume fraction is high, inertial effects are 
non-negligible and particles have finite size, i.e. size of the order of the relevant 
flow structures Q. 

In this study we focus on non-colloidal suspensions, mixtures where the dis¬ 
persed particles are greater than colloidal in size and thermal fluctuations are 
negligible. As Brownian motion is negligible there is no diffusion to create an 
equilibrium structure making the problem one of fundamental non-equilibrium 
physics. The aim of this study is to gain physical understanding of the role 
of the fluctuations induced by the suspended phase and their coupling to the 
mean flow, the effect of particle inertia and the modifications of the particle 
interactions when increasing the (bulk flow) Reynolds number. As shown also 
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here, it is fundamental to examine the local particle concentration, migration 
and segregation for a full comprehension of the transport processes at work. 
Inhomogeneities in the particle distribution are documented at low and finite 
Reynolds numbers, e.g. the so-called Segre-Silberberg effect Here we docu¬ 
ment how the interactions between the turbulent flow structures and particle- 
induced disturbances alter the macroscopic flow behaviour. 

Only few studies have been devoted to the inertial flow of suspensions in the 
presence of finite size particles. Matas et al. ^ performed experiments with a 
suspension of neutrally buoyant particles in pipe flow and defined the laminar 
and turbulent regimes according to the spectra of the pressure fluctuations be¬ 
tween the inlet and exit of the pipe. The critical Reynolds number separating 
the existence of the two regimes exhibits a non-monotonic behaviour with the 
volume fractions for large enough particles. A result partially reproduced by 
the numerical simulations in Yu et al. . Since velocity fluctuations exist at all 
Reynolds numbers, these authors choose the streamwise velocity perturbation 
kinetic energy as the criterion to distinguish between laminar and turbulent 
flow. A more detailed study on the transition of finite-size particle suspensions 


is performed by Loisel et. al. [6| for a fixed volume fraction of about 5%. The 
observed reduction of the critical Reynolds number is explained by the break¬ 
down of the coherent flow structures to smaller and more energetic eddies, which 
prevents the flow re-laminarization when decreasing the Reynolds number. The 
characteristics of a fully turbulent channel flow laden with finite-size particles 


are presented in [^, such as the decrease in the Von Karman constant with 


increasing volume fraction and the increase in the overall 


drag. 


on the inertial flow 


The present work extends the analysis of Lashgari et al. 
of suspensions of finite-size neutrally buoyant spherical particles. In this previ¬ 
ous study, we document the existence of three different regimes when varying 
Reynolds number. Re, and particle volume fraction, $. A laminar-like regime 
where viscous stress exhibits the strongest contribution to the total stress, a 
turbulent-like regime where the turbulent Reynolds stress mainly determines 


the momentum transfer across the channel (see also Q) and a third regime. 
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Figure 1: The instantaneous particle arrangement for (a) a turbulent-like flow, Re = 5000 & 
■F = 0.1, and (b) a particle-stress dominated flow, Re = 2500 &: <!> = 0.3. The data are not 
visualized at true scale. 


denoted as inertial shear-thickening, characterized by a significant enhancement 
of the wall shear stress that is not due to an increment of the Reynolds stress 
but due to the strong contribution of the particle stress. In the present work, 
we move our attention from the bulk flow behavior to the local behavior by 
studying in detail the particle dynamics, single and pair particle statistics. In 
particular, we examine the particle local volume fraction, dispersion coefficients 
and collision kernels for the three regimes introduced in Lashgari et. al. 

Our dataset is based on fully resolved numerical simulations of the two-phase 
system. 

We aim to connect our results to the seminal work by Bagnold B- Using 
experimental data of a suspension of neutrally buoyant solid particles in an an¬ 
nular domain between two concentric cylinders, Bagnold understood that the 
shearing of closely spaced particles would generate a normal or dispersive stress 


in addition to the shear stress [iB. He used the ratio between the grain inertia 
and the viscous stress to define different flow regimes. The viscous and iner- 
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tial regimes introduced by Bagnold are characterized by a linear and quadratic 
relation between the wall shear/normal stress and the shear rate, respectively. 
Inspired by Bagnold’s experiment, Fall et aZ. performed a similar study in 
plane Couette flow; these authors show a smooth transition from the Newtonian 
(viscous) to the Bagnoldian (inertial) regime by increasing the shear-rate. The 
laminar flow at high volume fractions behaves similarly to dry granular flows 
Q : the flow experiences discontinuous shear-thickening and fast particle migra¬ 
tion toward the regions of low shear. Both effects (shear-thickening and particle 
migration towards region with low shear) have been observed in several previ¬ 
ous investigations of dense suspensions at low Reynolds number, see |12l.ll3l.ll4| 
among others. Shear-thickening at higher volume fractions is examined among 
others in Haddadi and Morris |15j| who clearly identify the role of friction among 
particles in relative motion. The origin of shear-thickening in the presence of 
non-negligible inertia is attributed to the particle dynamic excluded volume in 


•H 


the recent work by Picano et a/.[16|l at lower volume fractions. The effective vol¬ 
ume fraction of the suspension increases because of the shadow region (a region 
with statistically vanishing relative particle flux) around the particles. Particle 
migration across the channel is not an inertial effect and is observed also in 
Stokes flow at high volume fractions [l^. The particles tend to migrate from 
regions of high to low shear due to the imbalance of the normal stress resulting 
from the particle interactions 171. 


Less is known of the inertial Bagnoldian regime. It is worth mentioning that, 
for the same bulk behavior, the Bagnoldian regime can be either Reynolds stress 
or particle stress dominated, as deduced from the data in Q. This finding moti¬ 
vated the present study where we focus on the particle dynamics to understand 
the two different underlying physical mechanisms. 

Understanding the dynamics of particle dispersion and collisions, especially 
when the particle inertia is non-negligible and the suspension is not dilute, is 
therefore important due to their direct connection to the flow bulk properties, 
as also demonstrated in this study. The mutual and hydrodynamic interactions 
between the particles produce irregular motions, promote lateral migration from 
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the instantaneous average particle trajectories and induce dispersion (for more 


details see 


la,llSj). As an example, we report in figured] the instantaneous 


particle distribution for two different regimes: i) a turbulent flow where trans¬ 
port is mainly determined by the Reynolds stresses and ii) a shear-thickening 
flow dominated by the particle stress. Note that the wall normal direction is 
amplified by a factor 5 for the sake of clarity and the particle colors represent the 
magnitude of their translational velocities. We note an uniform concentration 
for the turbulent-like flow (left panel) and an accumulation towards the channel 
centre for the flow dominated by the particle stress (right panel) that will be 
quantified and analyzed in this paper. 

Particle collisions are also relevant to the total momentum transfer and can 
be estimated from the relative position and velocity of the particle pairs H: 
these can be directly connected to the particle di ffus i vity in the cross-stream 
direction and to accumulation in specific regions 


21 


22 


23|. The opposite is true 


for Brownianan suspensions where the particle concentration variation arises 
from gradient-induced diffusivity [l^, and finite-size effects are less important. 
In this work we show a strong shear-induced self-diffusivity at high particle 
volume fractions which is not dependent on the Reynolds number and plays 
an important role in the collision dynamics and eventually on the bulk flow 
behavior. 

This paper is organized as follows. We discuss the governing equations, the 
numerical method and validations in §2. The results of the simulations are 
discussed in §3, whereas conclusions and final remarks are presented in §4. 


2. Governing equations and numerical method 

2.1. Governing equations 

We study the motion of suspended rigid neutrally buoyant particles in a 
Newtonian carrier fluid. The Navier-Stokes and continuity equations govern 
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( 1 ) 


the motion of the fluid phase, 

+ u • Vu) = —VP + + p-^f, 

V • u = 0, 


where p-^ and are the density and viscosity of the fluid. Here we denote the 
spanwise, streamwise and wall-normal coordinates as (x, y, z) with correspond¬ 
ing velocities u = {u,v,w). The force on the fluid, f, is due to the presence of 
the finite-size particles. The motion of the particles is governed by the Newton- 
Euler equations 


dt 


= F^ 


dflP 

TP ^^C _ rp,p 

dt 


( 2 ) 


where and 1^ are the mass and momentum inertia of particle p, and 
HP the velocity and rotation rate and and the net force and momentum 
resulting from hydrodynamic and particle-particle interactions that for neutrally 
buoyant particles read 


pP = f [-PI -H pf (Vu 4- Vu^)] • ndS + Fc, 


dV„ 


( 3 ) 


TP = f r X {[-PI + pf (Vu + Vu^)] • n}dS' -h Tc. 
JdV„ 


In this expression dVp represents the surface of the particles with unit normal 
vector n. The radial distance from the centre of the particle is denoted by r 
and the force and torque resulting from particle-particle (particle-wall) contacts 
are indicated by Fc and Tc. The equations for the fluid and particle phase 
are coupled by the no slip and no penetration conditions on each point X on 
the surface of a particle, i.e. u(X) = U^(X) = U(( -I- x r. We use here an 


Immersed Boundary Method 


24| . where this condition is satisfied indirectly by 


applying the forcing f on the right hand side of the Navier-Stokes equations. 


2.2. Numerical method 

We employ a Navier-Stokes solver coupled with an Immersed Boundary 
Method (IBM) to follow the motion of the fluid and rigid spheres in the domain. 
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The direct forcing method was originally proposed by Uhlmann 2^ and modi¬ 


fied by Breugem 


24j to ensure second-order spatial accuracy. An Eulerian fixed 


mesh is used for the fluid phase and a Lagrangian mesh to represent the moving 
surface of the particles. The IBM forcing imposes no-slip and no-penetration 
boundary conditions on the surface of the particles. When the distance between 
particles or with a wall becomes smaller or of the order of the mesh size, the 
interactions between the particles include an additional lubrication correction. 
Surface roughness effects are accounted for at very close approach. Finally, a 
soft-sphere collision model is employed to model collisions/contacts from the 
relative velocity and (slight) overlap of colliding particles, (see the appendix of 


26l| for more details). The IBM code has been used to study passive and active 


suspensions by Lambert et al. 


26| and Picano et al. [161, a. 


2.2.1. Validation 


The accuracy of the code has been verified against several test cases in 2^ . 
In particular, the benchmark cases include: 

I) Flow through a regular array of spheres at solid volume fraction of 6.5% for 
which 2"'* order accuracy has been shown together with excellent agreement 
with theoretical predictions in 271. 


2) Lubrication force between 2 spheres: accurate predictions up to a gap width 
of 2.5% of the sphere radius at resolution D/dx = 32 (see Q and accuracy test 


m- 


at 1 specific gap width in 

3) Drafting-kissing-tumbling interaction between two spheres for which 2"'^ or¬ 
der accuracy is confirmed. 

4) Validations of the global suspension behavior are reported in H3- 

5) The code has also been compared with the experimental data by Cate et. 
al. on a falling spherical particle in a closed rectangular container. The 
container dimensions were as in the experiment: 100 x 100 x 160 mm in the 2 
horizontal and the vertical direction, respectively. The diameter of the nylon 
sphere was 15 mm. The mass density of the nylon sphere 1120 kg/m3. The 
mass density of the working fluid, silicon oil in this case, was 960 kg/m3 (so the 
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Figure 2: Particle settling velocity versus observation time; the solid line represents the nu¬ 
merical results whereas the dots are the experimental data from 


mass density ratio = 1.17). The dynamic viscosity of the working fluid is 0.058 
kg/(ms) for estimated Reynolds number Re ~ 31.9. Figure [2] shows the particle 
velocity as a function of time; the dots are the experimental data. The small 
difference at early time is probably related to the different release procedure in 
simulations and experiments. 

In this work, we further validate the code by comparing the trajectory of 

a particle pair in homogeneous shear flow against the work of Kulkarni and 

Morris [^. These authors employ a Lattice-Boltzmann method to study the 

effect of particle inertia in a laminar flow. We choose the same box size of 

20a X 20a x 20a with a the particle radius and place the origin of the coordinate 

system at the center of the box. The two particles are initially positioned at 

Xi = (0,—4.85a, 0.32a) and X 2 = (0,+4.85a,—0.32a) and move in opposite 

• 2 

direction. The particle Reynolds number Re^ = = 0.1 with 7 the imposed 

shear rate and v the fluid kinematic viscosity. The initial particle velocity is 
the same as the local fluid velocity at the center of the particle and the initial 
rotation is equal to half the local fluid vorticity. The trajectory of the particle 
centers are displayed in figure [3] together with the data from [30j (open circles). 
The arrows display the direction of the particle velocity at the initial time. The 
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Figure 3: Trajectory of a particle pair released in a simple shear flow from initial positions 
Xi = (0,—4.85a, 0.32a) and X 2 = (0,+4.85a,—0.32a). The particle radius a = 0.05. The 
solid line indicates our simulation while open circles are the reference data in . The arrows 
show the initial direction of the particle velocities. 


results show a good agreement: the particle trajectories deviate away from the 
centerline after their interaction (a typical inertial effect). 


2.3. Flow configuration and numerical setup 

We simulate a channel flow with periodic boundary conditions in both stream- 
wise and spanwise directions. The box-size is 6h x 2h x 3h in the streamwise, 
wall-normal and spanwise directions. The particles have all the same radius, 
a = h/10. The Reynolds number is defined as Re = , where Ub and v 

are the fluid bulk velocity and kinematic viscosity. A wide range of parameters 
have been considered; the Reynolds number 500 < Re < 5000, and the particle 
volume fraction 0 < <i> < 0.3. The simulations start with a high amplitude 
localized disturbance in the form of two counter-rotating streamwise vorticities 
to efficiently trigger turbulence, if the Reynolds number is high enough for 
it to be sustained [32|. The particles are randomly positioned at time zero, with 
velocity equal to the local fluid velocity and rotation equal to half the local fluid 
vorticity. 
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For the simulations presented in this work we employ 480 x 160 x 240 uniform 
Eulerian grid points in the streamwise, wall-normal and spanwise directions and 
746 Lagrangian grid points to cover the surface of each particle. The statistics 
are computed when the flow is fully developed. To check if the results are statis¬ 
tically converged, we repeat the analysis using half of the number of the samples 
and compare the outcome with the one from the total number of samples. 


3. Results 

In this work we study the local properties of the particulate flow, in particular 
particle distributions, dispersions and collisions, and connect the results with 
the bulk flow regimes identified in [sj. 


3.1. Inertial regimes and Bagnold theory 

In this section, we analyse in detail the momentum budget of the two phase 


B- 


flow and provide a comparison with the seminal work by Bagnold l9|. Batch¬ 
elor 0 was probably the first to derive an analytical expression for the bulk 
stress of suspensions of rigid particles and to discuss the relation between the 
macroscopic properties of a homogenous suspension and the flow structures at 
the particle scales. He assumed that i) the bulk stress depends on the instan¬ 
taneous particles configuration in a flow element containing a large number of 
particles; ii) the configuration in each element depends on the history of the 
motion (memory effect). This shows the importance of the local microstructure 
and of its time history to determine the bulk behavior of the suspension. For 
colloidal suspensions in the inertialess regime, the relation between the particle 


and bulk scale structures are thoroughly reviewed by Morris 


0 . 


We employ the phase-averaged momentum equations following the formula¬ 
tion developed in |35l . l36l . l37 1 where the effect of spatial non-uniformity over a 
finite scale larger than the particle size has been taken into account, unlike in 
the original formulation by Bachelor. The phase average momentum equation 
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ids. 


(4) 


on the volume V with boundary 5(V) reads [37 1 

p ( + (1 - 0®^) dV = (f [^<7^ + (1 - 0®^] 

Jv Jsiv) 

where ^ is the phase indicator with values ^ = 0 and f = 1 for the fluid and 
particle phase. The stress and acceleration of the fluid and particle phase are 
denoted by and respectively. Assuming statistical homogeneity 

in the streamwise and spanwise directions, one can obtain an expression for 
the stress budget across the channel (see for the detail of the derivations the 
appendix of Q) 


dVf 


-[z/h) = -{v'*w'*) + + -{crL) = 


dVf 




(5) 


dz p ^ dz 

where T{z/h) is the total stress as a function of wall normal coordinate, z, nor¬ 
malized by the channel half width, h. The first term in the stress budget is the 
total Reynolds stress, = {v'*w'*) = (1 — consisting of 

the fluid and particle Reynolds stress weighted by local particle volume fraction, 
ip. The second therm, ry = y(l — is the viscous stress whereas the third 

term, Tp = —{(jP), is the stress due to the particles. The sum of the three 
terms is a linear function across the channel, as for a classic turbulent flow 


with the wall shear stress t„, = \ . The expression for the particle stress 

33|, the particle stress in the 


is discussed in detail in 


33, 


39 


37 1. Following 


absence of an external torque reads 


p - 1 V 
Y ^ 


[ + ajkXijnkdA - -^Ey [ ^-plf-xj + fjXijdV cr^-, 

Ja„ ^ y Jv„ ^ 


( 6 ) 


where x is the material point, /( is the local acceleration of the particle relative 
to the average acceleration, Ap and Vp are the particle surface area and volume 
and aij = —PSij + is the fluid stress tensor. The first term on the 

right hand side is the hydrodynamic stresslet resulting from the symmetric part 
of the shear stress. This is related to effective viscosity in dilute suspension |36|. 
The second term is the stress related to the particle acceleration and rotation 


12 















Figure 4: Profiles of (a) the Reynolds stress and (b) the particle stress across the channel for 
the cases indicated in the legend. 


with respect to the neighbouring flow and the last term is the inter-particle stress 
which depends on the near-held particle Interactions and particle collisions. 

The wall-normal prohles of Reynolds and particle stresses, tr and rp, are 
displayed in hgure 0] for hows at low and high particle volume fractions, $ = 
0.05 and <I> = 0.3. The cases at $ = 0.05 and Re > 2500 are characterized 
by a dominant contribution of the Reynolds stresses with slgnihcant particle 
stress only in the near-wall region. At high volume fractions, $ = 0.3, on the 
contrary, the particle stress accounts for more than 75% of the total stress and 
the Reynolds stress Is slgnihcant only in the intermediate region between the 
wall and centerline. As we will discuss in detail, the particle accumulation in 
the core region and in the layer close to the wall explains the high particle stress 
In the dense suspensions. Finally, we note that the major contribution to the 
total stress Is due to the viscous forces at $ = 0.05 and Re = 1000; Increasing 
Reynolds number at hxed $, the Reynolds stresses increases sharply when the 
how becomes turbulent while the particle stress slightly decreases. 

To understand the role of the different transport mechanisms on the bulk 
how behavior in the range of Reynolds numbers and particle volume fractions 
investigated, we show in hgure [5] (panels a,b,c) maps of the relative contribution 
of viscous, Reynolds and particle stress to the total momentum transfer inte¬ 
grated across the channel. The dashed lines represent Iso-levels of 25%, 50% 
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Figure 5: Contour map of the percentage contribution of (a) viscous stress, (b) Reynolds 
stress and (c) particle stress to the total momentum transport integrated across the channel. 
The isolines show the boundary of the regions in the map where the contribution of each term 
is more than 25%, 50% and 75%. (d) Contour map of effective viscosity, the normalized wall 
shear stress divided by the shear at the wall of the corresponding laminar flow. 
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and 75% of the total stress. The region where the viscous stress is more than 
50% of the total stress is limited to Re < 1900 and $ < 0.13. In this region, 
the action of viscous dissipation overcomes inertia. The contour lines in figure 
[5]Jb) show the non-monotonic behavior of the Reynolds stress which is also an 
indication of the level of fluctuations in the flow. This trend is in agreement 
with previous experimental and numerical findings P, [sj, where the authors re¬ 
port a non-monotonic behavior of the critical conditions for the occurrence of 
turbulent flow when increasing the particle volume fraction. The contribution 
of the Reynolds stress is more than 50% of the total for Re > 2000 and $ < 0.1: 
the fluid and particle phases induce strong fluctuations that cannot be damped 
by viscous dissipation. The region with $ > 0.13 is characterized by values of 
the particle stress larger than 50% of the total stress. In this region, we expect a 
high level of hydrodynamic and particle-particle interactions that induce strong 
particle stresses. Based onjSJb) andjSfc), the rate at which the particle stress 
contribution is increasing with Re is similar to the rate at which the Reynolds 
shear stress contribution is increasing (the lines have similar slopes). This sug¬ 
gests that for high $ the flow will not be dominated by turbulent transport 
when increasing the Reynolds number Q 

Finally, we display the effective viscosity in figure [SKd) . Following previous 
literature, we define the effective viscosity as the normalized wall shear stress 
divided by the shear at the wall of the corresponding laminar flow, /tq 

1^ . We observe a monotonic increase of the dissipation when increasing both 
the Reynolds number and the particle volume fraction. The regions where the 
contribution of each stress term is more than 50 % of the total (See panels 
a,b and c) are depicted on the map of panel (d) by solid blue lines. Following 
Lashgari et al. these regions represent the laminar, turbulent and inertial 
shear-thickening regimes where the viscous, Reynolds and particle stress con¬ 
tribute the most to the momentum transfer. The three transport mechanisms 
coexist with different relevance depending on the Reynolds number and particle 
volume fraction. 


To continue our analysis, we first recall the work by Bagnold [9j on inertial 
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Figure 6: Effective viscosity, the normalized wall shear stress divided by the shear at the wall 
of the corresponding laminar flow, fir = TwItq, vs Bagnold number for four representative 
values of the volume fraction (j). 


suspensions. This author introduces what is known as the Bagnold number, 
Ba = ARepy/X, where A = (q 74 /$)i/ 3 _i is the linear concentration computed 
as the ratio between the particle diameter and the average radial separation 
distance and Rep, the particle Reynolds number. This non-dimensional param¬ 
eter represents the ratio between inertial and viscous stresses and it is shown to 
describe the bulk behavior of the suspensions reasonably well. Bagnold defines 
a macro-viscous regime, Ba < 40, at low Reynolds number and low particle 
volume fraction where the relation between the stress and shear-rate is linear 
(similar to a Newtonian flow). The Bagnoldian regime, instead, Ba > 450, 
appears at higher Reynolds numbers and particle volume fractions and is char¬ 
acterized by a quadratic dependence of the stress on the shear-rate. 

The effective viscosity pertaining our simulations is depicted versus the Bag¬ 
nold number in figure [6] for four representative values of the particle volume 
fraction. Interestingly, we observe that the effective viscosity of the suspension 
is almost constant when Ba < 40, as expected in the visco-macro regime where 
the shear stress depends linearly on the Bagnold number (constant effective vis¬ 
cosity). When Ba > 70, all the curves collapse on a single line (see figure [5jl) 
and the effective viscosity varies linearly with the Bagnold number. 
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These results suggest that both the Reynolds and particle stress dominated 
flows fall into the Bagnoldian inertial regime: the same value of the Bagnold 
number and the same dissipation (effective viscosity) may therefore be explained 
by two different underlying physical mechanisms. The shear stress due to the 
residual turbulence becomes negligible when increasing the grain concentration, 
as predicted by Bagnold, and the particle stress takes the place of the Reynolds 
stress in the transport of momentum across the channel. At very high volume 
fraction the dynamics of the flow resemble granular media where the effect of 
interstitial fluid is negligible and the inter-particle collision is the main transport 
mechanism 41|. 


3.2. Single particle statistics 

The single particle statistics are computed by considering quantities related 
to each individual particle and taking ensemble average over time and space. In 
particular, we extract the local volume fraction, mean and rms velocities of the 
particles as a function of the wall-normal coordinate z. 

The wall-normal profiles of the local particle volume fraction are shown in 
figure [7] for different Reynolds numbers and particle volume fractions, covering 
the three different regimes introduced above. Based on the phase diagram in 
figure [HKd), we see that in the laminar regime the particles accumulate in the 
intermediate region between the wall and the channel centerline, 0.2 ^ z/h < 0.8 
( cf. data for Re = 500^ = 0.05 and d* = 0.1). This appears to be due to the 
Segre-Silberberg effect [3|, an inertial effect (inertial migration) resulting from 
particle-fluid interactions, in particular explained in the dilute regime by the 
balance between the Saffman lift 4^ and inhomogeneous shear rate and wall 
effects (see paper by Schonberg and Hinch The Segre-Silberberg effect is 

documented also at finite Reynolds numbers and volume fractions in the work by 


Matas et al. 


44| where it is shown that the particle equilibrium position moves 


closer to the wall when increasing the Reynolds number, i.e. increasing inertial 
effects. 

The data in the figure show that increasing the Reynolds number while keep- 
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Figure 7: Wall-normal profiles of the local volume fraction for different values of the Reynolds 
number Re and volume fraction see legend. 


ing a low volume fraction the flow becomes turbulent and the particle distribu¬ 
tion is almost uniform across the channel, except in the near wall region due to 
the one sided wall-particle interactions. The particle distribution is homogenized 
by the action of the Reynolds stresses for the flow cases at Re = {2500, 5000} 
and $ = {0.05,0.1}, which are in the turbulent region of the phase diagram in 
figure [5{d). 

The particle distribution in the inertial shear-thickening regime, i.e. $ = 0.3, 
exhibits a completely different behaviour with a significant accumulation of the 
particles in the core region. The tendency of the particles to migrate toward 
the channel centreline is not per se an inertial effect (shear-induced migration) 
and is attributed to the imbalance of the normal stresses in the wall-normal 
direction [see Ref. Q, [l^ for more details]. When particle layers are sheared 
over each other, normal particle stress tends to push the particle layers further 
apart from each other which causes migration towards the core. Indeed, Notte 


and Brady [45j | first documented particle migration towards the core region 
at Re = 0. This effect is also evident from the data presented here as the 
prohle of the local volume fraction does not change considerably increasing the 
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Reynolds number from 500 to 5000. The local particle volume fraction in the 
core region approaches the value for a random loose packing; here the particles 
experience an almost uniform translational velocity, of the same magnitude as 
the carrier flow in some sort of plug flow, see discussion below. We find a 
peak of the local particle volume fraction at zjh = 0.1 also at the highest 
volume fraction under investigation, corresponding to particle layering at the 
wall. Once a particle approaches the wall it tends to stay there because the 
interaction with neighbouring suspended particles is asymmetric and the strong 
near-wall lubrication force hinders departing motions. 

We hnally note that, as we increase the particle volume fraction above $ = 
0.2 at fixed Reynolds number. Re = 500, the particle distribution changes from 
that typical of the Segre-Silberberg effect (due to the fluid-particle interactions) 
to display a significant accumulation in the core region due to the particle- 
particle interactions (shear-induced migration). 

The mean and rms particle velocities are depicted in figure H] for the same 
representative cases above. Figure [Ua) shows the mean streamwise velocity 
component to highlight the slip velocities at the wall. The velocity profile is 
closer to the parabolic single-phase profile for the laminar cases, while it becomes 
blunt in the turbulent and inertial shear-thickening regimes. Increasing the 
Reynolds number, the mean flow becomes more uniform across the channel. 
The fluctuation (rms) velocities, figure [5Kb, c,d), have also nonzero values at 
the wall: the level of fluctuations is higher for the turbulent cases and lower 
for the laminar flows, as expected. Interestingly, the level of particle velocity 
fluctuations of the inertial shear thickening regime (in particular for the cases 
the cases Re = 2500, $ = 0.3 and Re = 5000, = 0.3) is similar to that of the 

turbulent flows closer to the wall, z/h < 0.35, and closer to the laminar values 
towards the centerline, zlh> 0.65. This observation suggests that the dynamics 
of the inertial shear-thickening flows are similar to that of the laminar flow in the 
core of the channel and to that of the turbulent flow in the near-wall region.The 
low level of fluctuations in the core region results from the significant particle 
accumulation discussed above, which also explains the large particle stress in 
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Figure 8: Wall-normal profiles of the mean particle streamwise velocity component and of the 
3 particle r.m.s. velocities for the different values of Re and $ indicated by the legend. 


the inertial shear-thickening regime. Finally, we note that the peak of the wall- 
normal velocity fluctuations is located at zjh = 0.1 for all the cases (except 
Re = 500 and $ = 0.05), indicating the relevance of the impact with the wall. 
The peak moves to zjh = 0.2 for the streamwise and spanwise components of 
the velocity fluctuations. 


3.3. Particle dispersion 

To understand the dynamics of the interactions in the flow, we study the 
particle dispersion. The hydrodynamic and particle-particle interactions in¬ 
duce lateral forces and diffuse the particles away form their initial path. Th e 
dispersion is quantified by the variance of the particle displacement [21 


Here, we compute the mean-square displacement of the particle trajectories, 
2 /p (^)7 a function of the time interval At, averaging 

over all times, t, and particles, p, sampling after the initial transient in the flow 
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development, 


< AXliAt) >=< [Xp(t + At) - Xp(t)]2 . (7) 


The sampling time for diffusion must be larger in the case of lower since the 
fewer interactions between the particles require more time to reach the statistical 
convergence [l^. The particle diffusion coefficients can then be calculated by 
measuring the slope of the mean square displacement 


D = 


< AXliAt) > 

2At 


( 8 ) 


where D is the matrix containing the six independent diffusion coefficients Dij. 
The normalized particle correlation is defined by 


R(At) 


< Xp(t + At)Xpit) >p,t ~ < Xp 
< X^ > - < Xp >2 


(9) 


and it is used to quantify memory effects. The particle correlation and mean 
square displacement can be directly connected by 


R(At) = 1 


-DAt 


( 10 ) 


The spanwise mean square particle displacement, defined by eq. o, is dis¬ 
played in figure m a) as a function of 7 At where 7 = ^ is the average shear-rate 
across the channel. For all the cases studied, the mean square dispersion grows 
initially quadratically. In this ballistic regime, the particle trajectories are cor¬ 
related and the displacements are proportional to At so that < Aa: >2oc (At)2. 
The trend changes for larger time lags when the classical diffusive behavior is 
retrieved. This is induced by particic-particle and hydrodynamic interactions 
that de-correlate the trajectories in time as shown among others by [l^ and 


46l|in Stokes and low Reynolds number flows. The asymptotic slope determines 


the dispersion coefficient: The highest mean square displacement are found for 
Re = 2500 — 5000 and $ = 0.05 — 0.1 in the turbulent regime. The strong fluctu¬ 
ations de-correlate the particles trajectories in shorter time and promote higher 
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dispersion. The lowest dispersion is obtained in the simulations at Re = 500 and 
$ = 0.05 — 0.1, i.e. in the laminar regime where the viscous stress is dominating 
the flow dynamics, see figure [5] 

The asymptotic trend of the mean square displacement of the cases at 
$ = 0.3 lies in between those of the turbulent and laminar regimes and smoothly 
increases with the Reynolds number. The two dashed black lines in the fig¬ 
ure represent therefore the lower and upper limit attained at Re = 500 and 
Re = 5000, the other cases lying in between. Thus, in the cases for which 
the particle stress provides the largest contribution to the momentum transfer 
(inertial shear-thickening), the diffusion coefficients are well below those of a 
turbulent flow. In the core region, the particle dispersion is low due to the sig¬ 
nificant particle accumulation, of the order of that in a laminar flow. However, 
the turbulent stress is still active (see figlHa) in the intermediate region between 
the wall and the centerline, which locally increases mixing. This explains why 
the overall dispersion in the inertial shear thickening regimes is between that of 
the laminar and turbulent regimes, as shown in figure |9l The Reynolds stress 
and particle stress dominated regimes exhibit therefore different values of the 
particle dispersion even if they assume the same value of the Bagnold number. 

The wall-normal mean square particle displacement is depicted in figure[9l(b). 
First, we note that the wall-normal dispersion is lower than its spanwise coun¬ 
terpart due to the presence of the walls. The data in figure IHKb) asymptotically 
approach the value 18^, the maximum possible wall-normal displacement when 
normalized with the particle radius. As observed for the spanwise dispersion, 
the highest and lowest diffusion pertain to the turbulent and laminar flows, 
whereas inertial shear-thickening flows display intermediate values. These data 
on single-particle dispersion provide additional evidence that two distinct dy¬ 
namics are at work in the turbulent or inertial shear-thickening regime, despite 
they can be both classified as inertial Bagnoldian flows. 

The values of the spanwise dispersion coefficients extrapolated for larger 
At from the data in hgure [9] are reported in table [TJ The largest dispersion 
coefficient is obtained for the case Re = 5000 and $ = 0.05 where the turbulence 
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Figure 9: Mean square particle displacement in the spanwise (a) and wall-normal direction 
(b) for the cases indicated in the legend. The thick black lines with arrows indicate the range 
of values assumed by flows in the inertial shear-thickening regime. 


23 




















Re 

500 

2500 

5000 

500 

2500 

5000 

500 

2500 

5000 

$ 

0.05 

0.05 

0.05 

0.1 

0.1 

0.1 

0.3 

0.3 

0.3 

Ox 

0.004 

0.095 

0.1 

0.005 

0.075 

0.085 

0.005 

0.018 

0.023 


Table 1: Dispersion coefficients computed by particle displacements in the spanwise direction. 


activity is the strongest and the lowest one for Re = 500 and $ = 0.05 in the 
laminar regime. As mentioned above, the wall-normal mean square displacement 
is limited by the walls and therefore its slope does not reach a constant value, 
corresponding to a well-defined value of the wall-normal dispersion coefficient. 


3.4- Particle-pair statistics 

We shall also consider the variations of quantities pertaining pairs of particles 
as a function of the distance between their centers, r. As r approaches the 
particle diameter, the near field interactions become important and collisions 
occur when r become less than one particle diameter. For the details of pair- 


particle statistics the reader is referred to the appendix of m among others. 
Here, we shortly introduce the Radial Distribution Function, g{r). In a reference 
frame with origin at the centre of one particle, the RDF is the averaged number 
of particle centers located in a shell of radius r and r -\- dr divided by the 
expected number of particles of a uniform distribution [see 22, Formally, 
g{r) is defined as 

1 dNr 


gir) = 


1 


47 r dr r^riQ ’ 


( 11 ) 


where is the number of particle pairs in a sphere of radius r and Uq = 
0.5 * Np{Np — 1)/V the density of particle pairs in the volume V, with Np the 
total number of particles. For small values of r, g(r) reveals the intensity of the 
particle clustering whereas g(r) —> 1 when r —>• oo (uniform distribution). 

The dynamics of the particle pair cannot be determined only by the pair 
distribution function. Following the study of Sundaram and Collins [^, we 
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compute the normal relative velocity of the particle pairs as function of r. Con¬ 
sidering particle i and j, the normal relative velocity of the particle pair is 
obtained as the inner product of their relative velocity and relative distance [see 
48| 


dVninj) = (Ui - Uj) 


Kr* - ^j)\ 


= (ui - Uj) • 


( 12 ) 


The normal relative velocity is a scalar quantity and can be either negative, 
dv~{rab) = dvnirab)\^Q, for approaching particles or positive, dv+{rab) = dwn(?'ab)| 
when the two particles depart from each other. The averaged normal rela¬ 
tive velocity can be therefore decomposed into < dvn{r) >=< dv^{r) > -|- < 
dv~{r) >. Finally, the collision kernel is obtained as the product of g{r) and 
< dv-{r) > 0 , 


K(r) = g{r)- < (r) > . 


(13) 


First we compute the Radial Distribution Function, g{r), the normalized 
probability of finding pair of particles at distance r. The results are shown in 
figure irol ai for several cases covering the three different regimes. For all the 
cases, the maxima of g{r) occurs at r/2a = 1 where the particles are in contact 
and the collision force is active. Note that to compute the RDF, r is discretized 
in the range of r/2a = [0.9,4] and the data for r/2a < 1 are not displayed (We 
may have some events where 0.9 < r/2a < 1 ). Increasing the distance between 
the particles, g(r) —>■ 1, as expected. The inertial shear-thickening cases are 
characterized by the highest value of g(r j2a = 1 ) and by an additional peak 
around r/ 2 a = 2 indicating the probability of finding a second layer of the 
particles. In all the cases, the peak of g{r) decreases slightly when increasing 
the Reynolds number as the inertia tends to decorrelate the particle paths. 

We then study the normal relative velocity as function of r, < du„(r) >. We 
show the statistics of the negative relative velocity, < dv~ >, in figure [TOTbi: 
this observable indicates the tendency of the particle pairs to approach each 
other. The relative velocity increases almost monotonically with r as the pairs 
are more likely to approach with higher speed when farther away. The highest 


> 0 ’ 
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values of < dv~ > occur for the laminar cases and the lowest in the inertial shear¬ 
thickening regime, also the densest cases (again the relative negative velocity 
increases smoothly with the Reynolds number for $ = 0.3). The turbulent 
and inertial shear thickening flows (the Bagnoldian inertial flows) both exhibit 
values of < dv~ > lower than the laminar flows. 

The physical mechanism for the reduction of < dv~ > is however different. 
While the turbulent flow tend to homogenize the suspension, the shear induced 
migration observed in the inertial shear-thickening regime produces a significant 
accumulation in the core region where the particles are transported downstream 
by the core flow at almost constant velocity. This suggests that inertia and 
turbulent eddies determine the local and bulk behaviour of the suspensions at 
low particle volume fraction, whereas particle interactions and shear-induced 
migration governs the behaviour of the flow at high particle volume fractions. 
The latter effect being almost independent of the flow inertia as discussed above. 

Using single particle statistics, we have shown that the particle dispersion 
is highest, moderate and low for the turbulent, inertial shear-thickening and 
laminar flows respectively. On the other hand, the analyses of the particle-pair 
statistics, i.e. relative velocity, shows the opposite ordering (laminar, turbulent 
and inertial shear-thickening regimes from high to low). This indicates that 
these two different aspects of the particle dynamics reflect the presence of the 
three different bulk regimes in a different fashion. 

The collision kernel, the product of g{r) and < dv~ >, is depicted in figure 
me)- Similarly to the normal relative velocities, this kernel increases monoton- 
ically with the distance r. Interestingly, the kernel assumes similar values when 
the particles are in contact (r/2a = 1) for all the cases studied. 

To explain the observation that the kernel function has similar values at 
contact, we display in figurefTTlthe separate contribution of g{r) and < dv~(r) > 
to the kernel at r/2a = 1, where the vertical axis are set to cover ±50% of the 
mean values of the data in each plot to ease a visual comparison. The values of 
g{r/2a = 1) and < dv~[r/2a = 1) > are almost independent of the Reynolds 
number; they do however increase and decrease with the particle volume fraction 
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Figure 10: (a) Pair distribution function, g{r)-, (b) Relative normal velocity, < dvn >; (c) 
Collision kernel K.(r) = g(r)- < dv^ (t) > as a function of the distance between the particle- 
pair, r, for different values of Re and 


in such a way that their product is almost constant and equal to 0.022± 10% for 
the collision kernel at r/2a = 1. This unique value seems to be valid for all the 
Reynolds numbers and particle volume fractions studied in the present work as 
well as those in [sj. Note that the error bar of these statistics is about 1%. 

As shown above, the dynamics of the particles changes considerably across 
the channel. To this end, we divide the channel into two regions: region /, close 
to the walls (0.05 < z/2h < 0.35 and 0.65 < z/2h < 0.95) and region II the 
middle of the channel (0.35 < z/2/i < 0.65). Note that the particle centers move 
in the range 0.05 < z/2h < 0.95, as the particle radius a = 0.05. 

We thus examine the radial distribution function, negative normal rela¬ 
tive velocity and the kernel operator in these regions, limiting the analysis to 
three flow cases representing the laminar, turbulent and inertial-shear thicken¬ 
ing regimes, see figure[T2l Note that the g{r) is normalized by the total number 
of particle pairs in each of the two regions. The normal relative velocity, see 
figure [T^ bl. is higher in the near-wall region for all the three cases, a fact at- 
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Figure 11: (a) Pair distribution function, g(r = 2a); (b) Relative normal velocity, < dvn > 
(r = 2a); (c) Collision kernel K(r = 2a) = g(r = 2a)- < dVn (r = 2a) > versus the Reynolds 
number for the indicated values of the volume fraction 


tributed to the strong background shear, however the difference decreases when 
the turbulent activity increases. The collision kernels, see fT^ cl. reveal that the 
kernel of the inertial shear-thickening flow, Re = 2500 and <5 = 0.3, is similar to 
that of the turbulent flow. Re = 5000 and <5 = 0.1, near the walls and to that 
of the laminar flow. Re = 500 and <5 = 0.05, in the flow bulk. This is inline 
with the results in figure [8] pertaining the particle velocity fluctuations. 

For the inertial shear-thickening flow, the radial distribution function at 
contact, g{rl2a = 1), is higher in region II. The difference between region I 
and II IS reduced for the turbulent flow, reflecting the more uniform particle 
distribution. The opposite behaviour is observed in the laminar regime where 
g(r j2a = 1) is instead slightly higher in region I. 

4. Conclusion and remarks 

We study the flow of suspensions of the finite-size neutrally buoyant particles 
in a channel, aiming to connect the local particle behaviour to the bulk flow 
properties. The analysis is based on data from direct numerical simulations 
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Figure 12: (a) Pair distribution function (b) Relative normal velocity (a) Collision kernel 
as a function of distance between the particle-pair, r, in the near-wall region I and channel 
core region II for three configurations. Laminar flow: Re = 500 & $ = 0.05, Inertial shear¬ 
thickening flow: Re = 2500 = 0.3 and turbulent flow: Re = 5000 & $ = 0.1 
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covering a wide range of Reynolds number, 500 < Re < 5000, and particle 
volume fraction, 0 < $ < 0.3, where the particles are rigid spheres with fixed 
ratio between the particle diameter and channel height of 1/10. 

The analyses of the stress budget reveals the existence of the three differ¬ 
ent flow regimes: laminar, turbulent and inertial shear-thickening depending 
on which of the stress terms, viscous, Reynolds or particle stress, is the ma¬ 
jor responsible for the momentum transfer across the channel. We show that 
both Reynolds and particle stress dominated flows fall into the Bagnoldian iner¬ 
tial regime Q: the suspension effective viscosity, i.e. the normalised wall shear 
stress, from the different simulations collapses when plotted versus the Bagnold 
number. Therefore, turbulent and inertial shear-thickening flows may share the 
same Bagnold number while the underlying momentum transport and dissipa¬ 
tion are distinct. 

Examining the particle distribution we show that in the viscosity dominated 
laminar flows, characterized by low particle volume fraction and Reynolds num¬ 
ber (<& < 0.1 and Re < 1000), the particles tend to accumulate at certain 
wall-normal equilibrium positions, a clear signature of the Segre-Silberberg ef¬ 
fect. The turbulent particle-laden flow, $ < 0.1 and Re > 1500, is instead 
characterised by a more uniform particle distribution due to the mixing by the 
turbulent eddies. At high volume fractions, $ > 0.2, we report a significant 
migration of the particles toward the channel centerline for all the Reynolds 
numbers under investigation, which explains the large contribution of the par¬ 
ticle stress in the so-called inertial shear-thickening regime. The particle accu¬ 
mulation in the core region is not necessarily an inertial affect as we observe a 
negligible variation of the local particle volume fractions when increasing the 
Reynolds number. 

The mean particle velocity profile becomes more blunt as the flow regime 
changes from laminar to either turbulent or particle-dominated shear-thickening. 
Interestingly, the velocity fluctuation amplitudes pertaining the inertial shear¬ 
thickening are closer to those of the turbulent flow in the near wall region, while 
they almost overlap to those of the laminar flow in the vicinity of the channel 
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centreline. The particle dynamics in the inertial shear-thickening regime appear 
therefore to share similarities with both the laminar and turbulent flow depend¬ 
ing on the wall-normal position. This is further conhrmed by examining the 
spanwise and wall-normal particle dispersion. For the inertial shear-thickening 
flows the turbulent activity is limited to the near-wall region while in the cen¬ 
tre of the channel the particles form a dense layer and are transported by the 
smooth flow of the carrier fluid. As a result of this, the dispersion coefficients 
of the inertial shear-thickening regime lie in between those of the two other 
regimes. Finally, we note that both the mean and fluctuating particle velocities 
exhibits slip velocities at the wall. 

We further consider the pair particle statistics in the three different regimes. 
In particular, we examine the pair distribution function, g(r), the approaching 
relative velocities, < dv~ > (r), and the collision kernel, K(r) = g(r) < dv~ >, 
as a function of the distance between the particle pairs. The laminar cases show 
the highest values of < dv~ > (r) and K(r) while the turbulent flows assume 
lower values due to the homogeneity created by the turbulent eddies. The 
lowest values of < dv~ > (r) and K(r) are found in the inertial shear-thickening 
regime, as a consequence of the particle packing in the core region. Separating 
the analysis in near-wall and centerline region, we observe larger values of < 
dv~ > (r) and K{r) in the wall region due to the strongest background shear. 

We have therefore demonstrated that the local particle dynamics clearly 
reflects the existence of the different flow regimes. The Bagnold number is 
shown to correctly predict the bulk flow behavior for the parameter range of 
this study; however, the details of the momentum transport and dissipation 
are different in the turbulent and particle-dominated regimes at low and high 
volume fraction, something which should be considered in any modeling effort. 
We believe future work should consider the role of the particle and fluid inertia 
and the particle dynamics in a mixture of particles of different size. 
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